Origin of Pyroelectricity in LiNbOa 
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We use molecular dynamics with a first-principles based shell model potential to study pyroelec- 
tricity in lithium niobate. We find that the primary pyroelectric effect is dominant, and pyroelec- 
tricity can be understood simply from the anharmonic change in crystal structure with temperature 
and the Born effective charges on the ions. This opens a new experimental route to studying py- 
roelectricity, as candidate pyroelectric materials can be studied with X-ray diffraction as a function 
of temperature in conjunction with theoretical effective charges. We also predict an appreciable 
pressure effect on pyroelectricity, which could be used to optimize materials pyroelectricity, and the 
converse electrocaloric effect, peak as T c is approached. 
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The theory of ferroelectricity had a classical period 
that culminated in the 1970s [l|-|3(, followed by a qui- 
escent period, and was rejuvenated in the 1990s with the 
introduction of modern electronic structure methods to 
these complex, interesting, and useful materials 0, HI- 
The fundamental physics of pyroelectricity, the change in 
polarization with respect to temperature, has not been 
re-investigated until now, and there is no previous first 
principles computation of pyroelectricity. Pyroelectricity 
is of current great interest since the discovery of particle 
acceleration of ions from changes in temperature at py- 
roelectric surfaces sufficient to generate hard X-rays in 
a commercial product as well as neutrons in heavy 
water via fusion Q . LiNb03 is a uniaxial pyroelectric 
with space group R3c in the polar phase with ten atoms 
per primitive cell, and a T c of 1480K 0, [loj]. The struc- 
ture, polarization and lattice dynamics of LiNb03 have 
been previously studied from first-principles using total 
energy, Berry's phase and linear response methods within 
Density Functional Theory (DFT) LiNb0 3 has 

been studied extensively experimentally [l|, [13 1 due to its 
use in SAW filters and non-linear optics. There is also 
much interest now in the converse of the pyroelectric ef- 
fect, the electrocaloric effect, for refrigeration or energy 
scavenging, [li- 16| 

Pyroelectricity is the change in spontaneous polariza- 
tion P s with temperature T. The total pyroelectric coef- 
ficient is 



„ dP s ,dP s . ,dP s ^ ,de s „ „ 

n = ^ = (^) e + (^M^). = n 1 + n 2 . (i) 



The first term on the right side is the primary pyroelec- 
tric effect and the second the secondary effect. Exper- 
imentally the pyroelectric effect is measured under the 
constraint of constant stress. The experimentally acces- 
sible or proper pyroelectric coefficient is due to the adia- 
batic current flow J due to a slow change in temperature, 
II' = where T is the change in temperature T with 
time t. The II' of an undamped sample can be expressed 



as 



n 



n' = iii + n 2 



n. 



(2) 



i measures the variation of spontaneous polariza- 
tion with respect to temperature at constant strain 
(clamped), which arises from changes in phonon occu- 
pations and anharmonicity. II2 is the result of crystal 
deformation where the strain caused by thermal expan- 
sion alters the polarization via the piezoelectric effect 
as n 2i = ajkCjkimd 



where the indices label coor- 
dinate directions 17j, repeated indices imply summa- 



tion, dn m are piezoelectric compliances, Cjkim are elas- 
tic moduli, and ay are the thermal expansion coeffi- 
cients. IT3 = 2a\P s is the difference from the total and 
proper pyroelectric coefficients [l|, where ol\ is the lin- 
ear thermal expansion coefficient of the plane perpen- 
dicular to the polar axis. II' can be measured with 
charge-integration or dynamic pyroelectric techniques [lj, 
whereas III , II2 , II3 cannot be measured directly. Un- 
derstanding the components of Iii , II2 , II 3 is crucial in 
studying pyroelectricity and its origin. 

We performed density functional theory (DFT) com- 
putations and fitted the results to a atomistic shell 
model. We used Density Functional Perturbation The- 
ory (DFPT) to compute phonons, effective charges, 
and dielectric constants. Wc performed first principles 
calculations with the ABINIT package (T9I ] within the lo- 
cal density approximation (LDA) [201 ] . Lithium Is and 
2s electrons, niobium 4s, 4p, 4<i, and 5s electrons, as 
well as oxygen 2 s and 2p electrons were considered as 
valence states. We constructed pseudopotentials using 
the OPIUM package [2l|. We used a kinetic energy cut- 
off of 45 Hartree and sampled the Brillouin zone using 
a 6 x 6 x 6 Monkhorst-Pack mesh of special k points. 
The results were carefully checked against previous all- 
electron [l(| and pseudopotential 11, 

The structural parameters of LiNbC>3 in 



computations, 
its polar 



ground state as functions of volume were obtained by 
relaxing the cell shape and atomic positions at seven vol- 
umes from 92.55 to 110.19 A 3 (-5 to 20 GPa from the 
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TABLE I. First principles calculation of structure, sponta- 
neous polarization P s , constant volume specific heat capacity 
C v , volumetric thermal expansivity a of LiNbOa using DFPT 
and the ABINIT code. 
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(10- 5 /K) 


DFT(OK) 
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0.86 






QHLD(300K) 


5.184 
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94.04 


3.59 


MD(300K) 
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0.63 




2.63 


Exp.(300K) 


5.151" 


13.876 a 


0.70-0.71" 


95.8 C 


3. 24-3. 83 d 



a [22|; b [2j,[23j; c [2j; d [25H27J 



TABLE II. First principles calculation of the elastic moduli c, 
piezoelectric strain constants d and piezoelectric stress con- 
stants e of LiNbQ 3 using DFPT and the ABINIT code. 





Smith et aZ. [27] 


Yamada [311 et al. 


Present 
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(xlO^N/m 2 ) 


(xlO^N/m*) 


(xlO^N/m^) 


Cll 


2.030 


2.03 


2.18 


C12 


0.573 


0.53 


0.68 


C13 


0.752 


0.75 


0.78 


Cl4 


0.085 


0.09 


0.15 


C33 


2.424 


2.45 


2.40 


C44 


0.595 


0.60 


0.55 


C66 


0.728 


0.75 


0.75 
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(xl0 -li C/N) 


(xl0~"C/N) 


(xl0 _11 C/N) 




6.92 


6.8 


8.12 


d22 


2.08 


2.1 


2.37 


d?,i 


-0.09 


-0.1 


-0.15 


ds3 


0.60 


0.6 


0.81 
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(C/m") 


(C/m") 


(C/m") 


eis 


3.76 


3.7 


3.722 


e 2 2 


2.43 


2.5 


2.317 


eai 


0.23 


0.2 


0.219 


e33 


1.33 


1.3 


1.718 



fitted equation of state). We also optimized the atomic 
coordinates for the paraelectric symmetry R3c with the 
same sets of optimized lattices. 

We computed the phonon frequencies using DFPT 
on a 4 x 4 x 4 grid of q-points at each of the seven 
volumes. The frequencies were interpolated onto a 
finer grid using short-range force constants [28j]. Quasi- 
harmonic Helmholtz free energies were obtained from 
these frequencies as functions of temperature and vol- 
ume. Isotherms were fitted to the Vinet equation of state 
& 

The polarization was computed for each of the seven 
volumes using the Berry's phase method with a 4 x 4 x 20 
mesh of k points. The results of the calculations were 
checked for convergence with respect to the number of 
k points and the plane wave cutoff energy. We obtained 
P s as the differences of the polarization between the po- 
lar and centrosymmetric LiNbC>3 at each volume, making 
sure we are on the same polarization versus mode coor- 
dinate curve 13011 . 
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FIG. 1. Molecular dynamics simulation results for LiNbOa in 
the NaT ensemble. The results for P — — 5, 0, 5 GPa are pre- 
sented by triangle, square, and cross line respectively. Exp.l is 
the experimental value of congruent sample in ref. @]. Exp. 2 
is stoichiometric sample in ref. Q|; Exp. 3 is experiment in 
ref. [32| . Note that our results from classical MD, so II' does 
not go to zero at zero temperature as required by quantum 
mechanics. 



The shell model approach has proven to be a computa- 
tionally efficient and confident methodology for the sim- 
ulation of ferroelectric perovskites, including bulk prop- 
erties of pure crystals, solid solutions and super lattices, 
and also surfaces and thin films properties 33j . In this 
model, each atom is represented by a massive core cou- 
pled to a massless shell, and the relative core-shell dis- 
placement describes the atomic polarization. The model 
contains 4th order core-shell couplings, long-range Ewald 
interactions and short-range interactions described by 
the Rydberg potential V(r) = (a + br)exp(r/p). The 
parameters were fit from the DFT and DFPT results of 
total energies, forces, stresses, phonon frequencies and 
eigenvectors, Born effective charges, and dielectric con- 
stants for a number of distorted and strained structures. 
We then performed classical molecular dynamic simula- 
tions with DL_POLY package 34 1. 

We computed the spontaneous polarizations P s during 
the MD simulations. NaT ensembles in MD simulations 
capture the evolutions of the system volume and shape 
corresponding to applied pressure, temperature. As a 
result, MD simulations in NaT ensemble allow us to com- 
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TABLE III. comparison of pyroelectric coefficient of LiNbC>3 
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FIG. 2. Pyroelectric coefficients. IIi (top- left), II2 (top-right), 
II3 (bottom- left) for P — —5,0,5 GPa, represented by trian- 
gle, square, cross line respectively. Bottom-right panel shows 
II2 directly calculated from the MD (square line) simulations 
and computed with formula (pentagon line). 



pute dP s /dT, the total pyroelectric coefficient II in Eq. 
[TJ We also performed MD simulations in the NeT ensem- 
ble (constant strain) and obtained ITi , and the difference 
gives II2. We computed II3 from MD NaT simulations. 

The MD simulations were carried out in a super cell 
with 8x8x8 primitive unit cells, which is 5120 atoms 
(5120 cores and 5120 shells). The NaT ensemble allows 
the shape and volume to change at constant stress, by 
which II and II' can be obtained. P s decreases with tem- 
perature and drops to zero at the phase transition to the 
paraelectric phase at 1200 K (Fig. [TJ) , which agrees well 
with the experimental value of 1430K [22] and 1480K Q. 

In order to understand the effects of volume error and 
the effects of pressure, we repeated the MD simulations 
and analysis at ±5 GPa. We found that P = 5 GPa 
reduces the volume about 3.6%, and P = —5 GPa in- 
creases the volume by 4.0%. T c is 1400 K and 1000 K for 
P = — 5 and 5 GPa respectively. 

We separately computed IIi from MD simulations in 
the NVT ensemble. The volume of the target tempera- 
ture T v was taken from the previous NaT simulations. 
MD simulations at T = T v , T„±10, T u ±20 K were carried 
out to calculate IIi at T v in Fig. [2] and Table Hill II i 
decreases with temperature and pressures, as does H2, 
calculated by II -III. We found n 2 = 13.5 /iC/m 2 K and 
n 3 = 17.6 /iC/m 2 K at 300 K. While lacking the direct 
and complete experimental data of all the coefficients of 
pyroelectricity, we estimate them as listed in in Table 
IIIII by combining the reported data of ref. 26, 27 1 and Q. 
There is good agreement between the experiment and 
present calculations. 

As a check, we calculated U2 in an alternative way 
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FIG. 3. The average value of z, u, v and w in MD. The P 3 
and II can be calculated using z, u, v, w and Born effective 
charge Z* (triangle line), compared with experiment (circle 
line) and direct MD results (diamond line). 



as n 2 = a jk c jk i m d 3 i m = 2e 3 iai + e 33 a 3 for LiNb0 3 , 
where e 3 i,e 33 are piezoelectric stress constants (Voigt 
notation), which are obtained by the hrst principles cal- 
culation at zero pressure and zero temperature as listed 
in Table HU Using ctj obtained from NaT simulations, 
we computed II2 (Fig. [5]) , which agrees with direct MD 
results at low temperatures up to 700K. 

We found that ITi is dominant and II 3 is small. The 
absolute value of both IIi and II2 increase rapidly with 
temperature as T c is approached. 

Up until now there has not been a clear exposition of 
the origin of pyroelectricity and the electrocaloric effect, 
but they are often considered as resulting from increas- 
ing polarization disorder with temperature. We find that 
the effects can be understood from the changes in crys- 
tal structure with temperature, as a simple anharmonic 
effect. We determined the average structural parameters 
z, u, v and w [IH from the average atomic positions in 
the MD simulations (Fig. [3]) . We computed the P s ver- 
sus temperature using these average positions with the 
Born effective charges Z* obtained from the DFPT com- 
putations, and P s — j| Z*Ti where is the «th ionic 
displacement along the polar axis from the centrosym- 
metric to polar structures, e the elementary charge and 
f2 the unit cell volume (Fig. [3]). The results show that the 
pyroelectric effect can be entirely understood in the clas- 
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sical regime above room temperature from the change in 
average structure with temperature, peaking at T c . Thus 
the anharmonic internal atomic rearrangement with re- 
spect to the temperature contribute the dominant part of 
the pyroelectricity. We find that pyroelectric coefficients 
could easily be obtained experimentally without electri- 
cal measurements, simply by studying changes in crystal 
structure with temperature, along with first-principles 
theoretical effective charges Z* . Good pyroelectrics and 
electrocaloric materials should have T c slightly higher 
than the operating temperatures. 
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